Whole genome demographic models indicate divergent effective population size histories shape contemporary genetic diversity gradients in a montane bumble bee

Abstract Understanding historical range shifts and population size variation provides an important context for interpreting contemporary genetic diversity. Methods to predict changes in species distributions and model changes in effective population size (N e) using whole genomes make it feasible to examine how temporal dynamics influence diversity across populations. We investigate N e variation and climate‐associated range shifts to examine the origins of a previously observed latitudinal heterozygosity gradient in the bumble bee Bombus vancouverensis Cresson (Hymenoptera: Apidae: Bombus Latreille) in western North America. We analyze whole genomes from a latitude‐elevation cline using sequentially Markovian coalescent models of N e through time to test whether relatively low diversity in southern high‐elevation populations is a result of long‐term differences in N e. We use Maxent models of the species range over the last 130,000 years to evaluate range shifts and stability. N e fluctuates with climate across populations, but more genetically diverse northern populations have maintained greater N e over the late Pleistocene and experienced larger expansions with climatically favorable time periods. Northern populations also experienced larger bottlenecks during the last glacial period, which matched the loss of range area near these sites; however, bottlenecks were not sufficient to erode diversity maintained during periods of large N e. A genome sampled from an island population indicated a severe postglacial bottleneck, indicating that large recent postglacial declines are detectable if they have occurred. Genetic diversity was not related to niche stability or glacial‐period bottleneck size. Instead, spatial expansions and increased connectivity during favorable climates likely maintain diversity in the north while restriction to high elevations maintains relatively low diversity despite greater stability in southern regions. Results suggest genetic diversity gradients reflect long‐term differences in N e dynamics and also emphasize the unique effects of isolation on insular habitats for bumble bees. Patterns are discussed in the context of conservation under climate change.


ent in the bumble bee Bombus vancouverensis Cresson (Hymenoptera: Apidae: Bombus
Latreille) in western North America. We analyze whole genomes from a latitudeelevation cline using sequentially Markovian coalescent models of N e through time to test whether relatively low diversity in southern high-elevation populations is a result of long-term differences in N e . We use Maxent models of the species range over the last 130,000 years to evaluate range shifts and stability. N e fluctuates with climate across populations, but more genetically diverse northern populations have maintained greater N e over the late Pleistocene and experienced larger expansions with climatically favorable time periods. Northern populations also experienced larger bottlenecks during the last glacial period, which matched the loss of range area near these sites; however, bottlenecks were not sufficient to erode diversity maintained during periods of large N e . A genome sampled from an island population indicated a severe postglacial bottleneck, indicating that large recent postglacial declines are detectable if they have occurred. Genetic diversity was not related to niche stability or glacial-period bottleneck size. Instead, spatial expansions and increased connectivity during favorable climates likely maintain diversity in the north while restriction to high elevations maintains relatively low diversity despite greater stability in southern regions. Results suggest genetic diversity gradients reflect long-term differences in N e dynamics and also emphasize the unique effects of isolation on insular habitats for bumble bees. Patterns are discussed in the context of conservation under climate change.

| INTRODUC TI ON
Populations in complex landscapes exhibit differences in diversity and connectivity driven by the distribution of habitable areas (Holderegger & Wagner, 2008;McRae & Beier, 2007;Storfer et al., 2007). Changes in genetic diversity from human-associated habitat alteration are an increasing concern, and factors like landscape modification and climate change have been associated with reduced genetic diversity and connectivity (Epps et al., 2006;Jha, 2015;Rubidge et al., 2012;Struebig et al., 2011). However, genetic variation is also influenced by the legacy of historical events that have driven fluctuations in effective population size (N e ). The impact of Quaternary climate variation on genetic diversity is well-recognized (Hewitt, 2000), and recent studies have more explicitly examined links between species range stability and genetic diversity (Carnaval et al., 2009;Garrick et al., 2021;Koch et al., 2018). Revealing how landscape changes drive declines and recovery in N e has thus become important for reconstructing evolutionary histories, but also for interpreting contemporary genetic diversity and making predictions about vulnerability to environmental change.
Genomes contain substantial information about historical demography. The development of methods to reconstruct historical demography using whole genomes even from single individuals (e.g., Li & Durbin, 2011;Schiffels & Durbin, 2014) has thus facilitated highresolution temporal analyses in wild populations (Mather et al., 2020;Nadachowska-Brzyska et al., 2016;Spence et al., 2018). Sequentially Markovian coalescent (SMC) methods (Mather et al., 2020), in particular, have gained favor for rapidly modeling N e trajectories because these methods can reveal demographic history for a population using only a single diploid genome sequence, although can also be applied with population sampling (Foote et al., 2021;Nadachowska-Brzyska et al., 2016;Taylor et al., 2021). Essentially, SMC models estimate the time to the most recent common ancestor for local genealogies at each locus in the genome as defined by ancestral recombination events (Li & Durbin, 2011;Schiffels & Wang, 2020). By providing estimates of N e changes over time, as opposed to providing a singlepoint estimate of diversity (e.g., current heterozygosity), SMC models are a powerful tool for examining how climate fluctuations may influence genetic variation (Taylor et al., 2021).
Thus, examining how populations historically responded to climate fluctuations would be useful for better-understanding threats from environmental change even for species that remain common (Dellicour et al., 2017;Koch et al., 2018).
Bombus vancouverensis is common in western North America Stephen, 1957) and has been a focus of several population genetic studies (Christmas et al., 2021;Jackson et al., 2018Jackson et al., , 2020Lozier et al., 2013Lozier et al., , 2016Lozier et al., , 2021. The population genetic structure of B. vancouverensis is well understood from prior analyses of mitochondrial sequences, microsatellites, and RAD-tag sequences Jackson et al., 2018;Lozier et al., 2011Lozier et al., , 2013. This study focuses on a latitudinal transect across westernmost B. vancouverensis populations (California, Oregon, Washington, USA; Figure 1), which provides an elevational gradient from zero to >2900 m above sea level and represents a morphologically and partly genetically distinct lineage . This region is interesting because these populations are weakly structured (global F ST = ~0.02) but exhibit a clear latitude-altitude gradient in genetic diversity, with diversity lower in southern high-elevation populations while northern populations appear more well-connected with higher heterozygosity (Jackson et al., 2018). Another interesting aspect of B. vancouverensis in the region is the presence of populations on islands in the Pacific Northwest. Islands are well known to accelerate genetic signatures of isolation in bumble bees (Goulson et al., 2011;Jha, 2015;Lozier et al., 2011) and can serve as a reference to test whether SMC models can detect effects of recent discrete fragmentation on genome variation.
In this study, we take advantage of existing whole genome data from prior studies Heraghty et al., 2020Heraghty et al., , 2022 to address new hypotheses relating to the origins of latitudinal genetic diversity gradients in B. vancouverensis. Grounded in prior knowledge about latitudinal patterns of genetic diversity and K E Y W O R D S bumble bees, effective population size, elevation, latitude, MSMC2, PSMC, sequentially Markovian coalescent, species distribution models

T A X O N O M Y C L A S S I F I C A T I O N
Biogeography, Conservation genetics, Ecological genetics, Population genetics structure across the region, we combine demographic models from representative whole genome sequences with temporal species distribution models to examine how the history of climatic fluctuations has shaped range dynamics and genetic diversity. We refine estimates of spatial variation in heterozygosity using the resolution afforded by whole genome data and then employ SMC modeling to estimate changes in N e over time using individual and populationlevel genome sequences. We also include a genome from an insular population to evaluate signatures of recent isolation in SMC models.
We use extensive occurrence data to model the B. vancouverensis range over the last ~130,000 years and examine how N e changes are associated with fluctuations in range size and stability. We investigate several related hypotheses.
(1) If current genetic diversity differences observed among populations can be explained by long-term Quaternary N e differences, N e in less diverse southern populations will be persistently smaller than in northern populations through time.
(2) If more recent N e declines can be detected with SMC models, insular genome data should demonstrate strong contemporary 2 | ME THODS

| Whole genome data
Whole genomes for B. vancouverensis were obtained from previously collected samples that were used for investigations of taxonomy, assembly of reference genomes, or analyses of environment-genome associations Heraghty et al., 2020Heraghty et al., , 2022 ( Figure 1, see Appendix S1, Tables S1 and S2). We combine these previously published whole genome data to conduct novel statistical analyses relating to the history of previously reported latitudealtitude genetic diversity clines in this species (Jackson et al., 2018;Lozier et al., 2011). All bees included here are female diploids and should represent distinct colonies (Jackson et al., 2018).  Table S2). Based on initial results using individuals (see Results), for some analyses mainland samples were grouped into three regional populations ( Figures 1   and 2, Tables S1 and S2; southern = southern California, midrange = northern California + southern Oregon, northern = northern Oregon + Washington).
We used the SFS for all samples in each population group for diversity calculations, but for plotting single population SFS, we also performed calculations for nine randomly selected individuals per region so that the SFS would have the same bin sizes. We determined the two estimates of θ (assuming θ = 3N e μ for haplodiploid bees, which have 1.5 allele copies per individual; Goulson et al., 2008) based on pairwise differences and segregating sites (θ π and θ W , respectively) and Tajima's D for each scaffold (Korneliussen et al., 2013). Values were scaled by the number of sequenced bases to obtain per site values. We compared diversity among southern, mid-range, and northern regions using linear mixed models conducted with r v4.0.5 (R Core Team, 2021) package 'lme4' 1.1-27 (Bates et al., 2015) with the region as a fixed effect and scaffold specified as a random effect. P-values were determined with 'lmerTest' 3.1-3 using the Swatterwhaithe method (Kuznetsova et al., 2017), and a likelihood ratio test was used to compare the region model against a random effect-only model. Individual heterozygosites were determined using a folded SFS (-doSFS and re-AlsFs commands, as above) from each sample (the "1" count bin in the SFS divided by total counts). For each region pair, we also estimated the folded 2D SFS and weighted F ST (not performed for San Juan Island genome with only one sample). Visualization of population structure was performed using PCA with pcANgsd v1.03 (Meisner & Albrechtsen, 2018).

| PSMC
The Pairwise Sequentially Markovian Coalescent (PSMC) model (Li & Durbin, 2011) is a widely used (Foote et al., 2021;Mather et al., 2020;Taylor et al., 2021) SMC methods that estimate the distribution of coalescent times between two unphased chromosomes from a single genome to infer N e over a series of intervals back in time. PSMC inputs were generated using sAmtools mpileup, with minimum mapping quality set to q = 20, base quality set to Q = 20, the mapping quality downgrading parameter set to C = 50, and the 500 kb scaffold bed file to specify target regions. Conversion to the consensus diploid genome was performed using bcFtools cAll, followed by conversion to fastq using the sAmtools 'vcfutils.pl vcf2fq' script, setting sites with minimum root-mean-squared mapping quality (Q) < 20, within 10 bp of an insertion-deletion, and below or above coverage thresholds (d or D, respectively) as missing data. Samples were sequenced to similar coverage so we specified a minimum depth of d = 10 and a maximum depth of D = 50 (see below for details on high coverage samples that required alteration of D).
Default settings for the atomic time interval parameter −p = "4 + 25*2 + 4 + 6" were designed for human analyses (Li & Durbin, 2011) but are expected to work well for many organisms (Patil & Vijay, 2021). We modified other parameters to ensure at least 20 recombination events per interval. We set initial -t (maximum TMRCA) at 7 and initial -r (ratio of scaled mutation to recombination rate θ/ρ) at 2, although the shape of outputs varied little with changes to these parameters, as expected (Mather et al., 2020). For plotting, we used a generation time of one generation per year and a direct mutation rate estimate of 3.6 × 10 −9 per site per year (Liu et al., 2017). Values converted to "real world" time and N e were extracted using the psmc_plot.pl script (modified to rescale using 3N e u) with the -R option and plotted using geom_step of 'ggplot2' 3.3.5 (Wickham, 2016) in r.
We also had the three higher coverage genomes (Tables S1 and

| MSMC2
Although we primarily present results from psmc, ensuring that results are comparable across methods can strengthen inferences (Mather et al., 2020). msmc2 is similar to psmc but can use population-level data (Schiffels & Wang, 2020). We focus on comparisons of N e for the mainland regions, as San Juan Island did not have population-level sequencing. Given similar individual trajectories for spatially proximate genomes in PSMC, genomes from some nearby sites were combined for MSMC2 to increase population samples (6-10 haplotypes each). Data files were prepared with minor modifications of scripts from J Rick (https://github.com/jessi caric k/msmc2_scripts). MSMC2 also required the mappability mask.
Separate VCF files for each scaffold and individual were created using sAmtools and bcFtools as above and excluding sites with <10x coverage. Mask BED files representing sequenced sites for each sample and scaffold required for MSMC2 were created using 'bamcaller.
py' in the MSMC tools repository. Input files for populations were prepared using VCFs, individual masks, and the masked reference genome using "generate_multihetsep.py" in MSMC tools. MSMC2 was run setting the -I parameter to only evaluate coalescence between haplotypes within individuals for each population (Schiffels & Wang, 2020) and using the default time interval parameter (−p = "1*2 + 25*1 + 1*2 + 1*3"); however, initial results suggested that the default -p time interval setting produced massive variance among populations and several unrealistically large N e estimates at the most recent time interval (<2500 years before present, or ybp). Such spikes are a typical sign of model overfitting (Schiffels & Wang, 2020) and we used a less complex setting of −p = "1*2 + 16*1 + 1*2" recommended by Schiffels and Wang (2020) to eliminate the N e spikes without notably altering curve shapes within the main period of interest. We converted outputs to real world time and N e for plotting as above.

| Distribution modeling
To predict changes in habitable areas during the Pleistocene and Holocene we used species distribution models from the presenceonly "maximum entropy" method mAxeNt (Elith et al., 2011;Phillips et al., 2006  . We removed localities likely to represent B. bifarius based on expected distributions  and by limiting the spatial extent to Longitude −150 to −110 decimal degrees and Latitude 35 to 65 decimal degrees. We only retained observations that had specific descriptive locality information for the record (e.g., at a finer resolution than county) to minimize georeferencing inaccuracy.
We obtained current and historical climate data (BIOCLIM variables, e.g., Hijmans et al., 2005) from 'rpaleoclim' 0.9 (Brown et al., 2018). We focused on variables for current, last glacial maximum (LGM, ~20,000 ybp), and last interglacial (LIG, ~130,000 ybp) time periods (Karger et al., 2017(Karger et al., , 2020(Karger et al., , 2021Otto-Bliesner et al., 2006). Given the scale of the study, we used a relatively coarse resolution of 10 arc minutes. We reduced highly correlated variables using hierarchical clustering (hclust in R) and selecting one variable We used the mAxeNt version implemented in the R package 'maxnet' 0.14 (Phillips, 2017). We used 'SDMtune' 1.1.5 (Vignali et al., 2020) to train and refine models. We first filtered the data to reduce spatial bias by removing points within two raster cells of each other (Kiedrzyński et al., 2017). This resulted in a filtered data set of 483 presence points. We randomly split these into 80% training (N = 386) and 20% testing (N = 97) points. Although the data remain more abundant in the contiguous USA and southern Canada, filtered occurrence records were much more even and largely represent the known distribution of the species (e.g., Stephen, 1957;Thorp et al., 1983). We randomly sampled 10,000 points as background observations and extracted predictor values at each presence and background point using 'raster' 3.5-2 (Hijmans, 2020). We trained initial models using the 'SDMtune' train function. We used the 'SDMtune' optimize Model function to generate models from training data with different combinations of regularization multiplier values (0.5-10) and linear, quadradic, product, and hinge feature classes, and compared the resulting models using AICc. The best model included a regularization multiplier of 1.0 and linear and hinge feature classes.
We trained models using selected parameters, plotted the receiver operating characteristic (ROC) curves, and calculated the area under the ROC curve (AUC) for test and training data sets. As a final evaluation of model performance, we used k-fold cross-validation using randomly generated sample folds (k = 5) and folds generated using "checkerboard" spatial partitioning (checkerboard1 function, aggregation factor = 10) (Radosavljevic & Anderson, 2014) from the 'ENMeval' 2.0.1 package (Kass et al., 2021). Cross-validation models were trained using the 'maxnet' algorithm (called with 'SDMtune') and evaluated with AUC (test AUC calculated using mean test AUC across folds). Final models were constructed using the selected tuning options and were spatially and temporally projected onto the current, LGM, and LIG predictor rasters with the default mAxeNt clog-log 0-to-1 scaling.
To evaluate how heterozygosity and the PSMC-inferred changes in population size may relate to changes in habitat suitability between current and LGM landscapes, we first calculated the "bottleneck size" for each genome as the difference between maximum preglacial N e and minimum glacial period N e from PSMC (see Figure S1). Next, we calculated a metric of "climatic stability" as the difference in predicted suitability between current and LGM Maxent rasters (positive values indicate less suitability at the LGM, negative values indicate greater suitability at LGM) and extracted values for each locality. Given the similarity of current and LIG maps, we only calculated the current-LGM difference; changes in suitability at genotyped localities between LGM and today can be interpreted as sites that fluctuate in abiotic suitability given the relative symmetry for LIG-LGM. As a second metric of stability, we evaluated the proximity of sample sites from areas that would have been suitable during the LGM (the "minimum distance to suitable LGM point", Figure S1).
We determined this distance using the distGeo function in the 'geosphere' 1.5-14 R package (Hijmans, 2019) (default WGS84 ellipsoid) to calculate the minimum geographic distance from each sample site to an LGM model raster cell with a predicted suitability ≥ the mean suitability of sites in the current species distribution model (≥ 0.882).
Distances were log-transformed.
We tested the effects of change in suitability at each sample site or distance to the nearest suitable LGM point on both heterozygosity and PSMC-bottleneck size using lme4 linear mixed effects models, specifying locality as a random effect. The San Juan Island genome was excluded because the population decline persisted and became most severe after the LGM (see Results) and determining the N e minimum attributed to the glacial maximum vs changes following glaciation was challenging, as was determining how to incorporate the effects of discrete isolation by water. Inferences for the high coverage genomes from northern and southern regions were essentially the same as moderate coverage samples (Figure 1c), suggesting that differences in PSMC trajectories were robust to sequencing depth. There was no overlap in bootstraps among samples except for some convergence near the LGM, again suggesting that differences inferred in the moderate coverage genomes were robust. Despite its northern provenance, the high coverage sample from San Juan Island had a historical trajectory similar to the southern genomes, except for the most recent interval that indicated a severe and well-supported recent bottleneck after 12,000-13,000 ybp (Figure 1c). Inferences were similar when variants were determined using the mappability-masked reference genome ( Figure S2), indicating that mappability was not a major factor.

| MSMC2
Major inferences about regional variation in long-term population sizes were generally supported by MSMC2 ( Figure S3). The historical expansion, a late Pleistocene bottleneck, and some degree of postglacial recovery paralleled PSMC. Likewise, the expansion peak N e and late Pleistocene bottleneck size were larger in northern populations, with the magnitudes of both declining toward the south.
Population N e trajectories generally matched the individual-genome PSMC plots, both for estimates of N e and the timing of size changes ( Figure S3).

| Parallels between summary statistics and SMC models
Demographic modeling was consistent with summary statistics from ANGSD in terms of diversity and possible nonequilibrium dynamics. Despite differences in the magnitude of N e fluctuations, the mean N e over time was greatest in the north and lowest in the south, and mean N e was highly correlated with individual heterozygosity, as would be expected given the theoretical relationship between N e and heterozygosity (Pearson's r = 0.96, p < .001; Figure 3a). The downsampled SFS were consistent with the signatures of expansion from PSMC inferences ( Figure S4).
The number of variants estimated decreased from north to south, with approximately 2.7, 2.5, and 2.3 M variant sites estimated in northern, mid-range, and southern regions, respectively, but F I G U R E 3 Correlation between mean N e over PSMC time intervals and individual heterozygosity for all moderate coverage mainland genomes for B. vancouverensis in CA, OR, and WA, USA (a). Relationships between niche stability at each genome coordinate (difference between current maxent suitability and LGM maxent suitability) and heterozygosity (b), niche stability and the size of the PSMC-inferred glacial period bottleneck (c), and the minimum distance from each sampling coordinate to a pixel with suitable LGM values (suitability values equivalent to the mean current suitability of genome samples shown) and the size of the PSMC-inferred glacial bottleneck (d). See Figure S1 for a schematic illustrating some of these metrics.

| Distribution models
The final Maxent model had a training and test AUC = 0.86 ( Figure S5), and both random and checkerboard spatial crossvalidation also produced AUCs in the 0.85-0.86 range for training and test data. Models from spatially-filtered observations performed suitably even with some sampling bias toward the southern part of the species range, and the map corresponds with prior models and range maps Thorp et al., 1983). The three models (current, LGM, LIG, Figure 1a In contrast to expectations for a stability-diversity hypothesis, heterozygosity was inversely related to the difference in raster suitability between time points (Figure 3b), suggesting that less diverse sites are more climatically stable, although this was not significant (Table S3). There was no significant relationship between the difference in suitability at each locality with the PSMC-bottleneck size during the last glaciation (Figure 3c; Table S3a). The clearest relationship detected was for the interglacial-to-LGM bottleneck size with the minimum distance to a suitable LGM point (Figure 3d, Table S3b): Larger Pleistocene bottlenecks were evident in northern genomes that are more distant from areas with average suitability during the last glacial period. The distance results remained significant when the mean Maxent suitability for all SDM presence coordinates (0.609) was used as a more conservative threshold for minimum distance calculations, rather than the more stringent threshold from genomic sampling sites (0.882) ( Table S3c). Given that heterozygosity is greatest in the north and these localities are most distant from suitable LGM raster cells, current heterozygosity was not associated with the distance-based measure of climatic range stability.

| DISCUSS ION
Combining temporal coalescent analyses from modest genome resequencing with paleo-SDMs proved useful for understanding variation in contemporary genetic diversity in B. vancouverensis.
Our objective was to better understand north-to-south gradients in genetic diversity throughout the Sierra-Cascade Mountain region of the western USA and to test whether such differences may be explained by long-term evolutionary processes or might require a more recent demographic explanation. Analyses suggest that all populations experienced fluctuations in N e alongside major range shifts. However, variation in genetic diversity can be explained by divergent population histories impacting the degree of these fluctuations, with more diverse contemporary populations exhibiting signatures of elevated N e throughout the late Pleistocene, particularly during periods of expansion when the species range was more continuous.
SMC models suggest that B. vancouverensis populations fluctuate at the same temporal scales throughout the latitude-altitude cline examined in this study. It is unsurprising that Quaternary climate dynamics drive concurrent demographic and genetic dynamics (Hewitt, 2000;Shafer et al., 2010), but it is interesting that the magnitudes of fluctuations differ across populations corresponding to the distribution of contemporary genetic variation. Using temporal species distribution models, we examined possible explanations for patterns of diversity and the magnitude of population size fluctuations (Carnaval et al., 2009;Koch et al., 2018). Studies that have investigated the relationship between Quaternary niche stability and genetic diversity have often found that regions with high stability harbor the retention of genetic lineages over time (Carnaval et al., 2009;Vasconcellos et al., 2019), but in other species, including other North American Bombus, the inverse relationship has been observed . Most sites in the B. vancouverensis range experienced a decline in habitat suitability at the LGM compared with today and the LIG. But in contrast to expectations under a hypothesis of greater stability driving retention of genetic variation, high-diversity northern populations experienced the largest loss in LGM suitability. Also somewhat counterintuitively given patterns of contemporary variation, northern genomes experienced much more substantial N e reductions through the last glaciation than less diverse southern genomes. Bottleneck magnitudes were larger in the northern populations that are more distant to locations that would have been suitable in the LGM, indicating that these populations would have had to disperse farther to track shifting climates. Although the distance to the LGM location with mean suitability of contemporary sites is not a perfect metric, as sites with below-average suitability could still have been occupied during the LGM, the strong relationship with PSMC-bottleneck size suggests it is a reasonable proxy to capture the overall greater habitat loss leading to a larger bottleneck in northern populations during this time.
We hypothesize that the counterintuitive inverse results relating to genetic diversity, historical bottleneck size, and niche stability are driven by the much greater N e of northern populations during climatic periods when they expand to become larger and more wellconnected to other parts of the range (e.g., LIG and today), while southern populations remain relatively isolated at high elevations even during favorable periods. Although most populations occupy areas that would have been unsuitable during the LGM, southern populations would have had less genetic variation to lose at this time owing to the weaker interglacial expansion. Furthermore, the southern B. vancouverensis range did not suffer a regional collapse of suitable habitats at the LGM; high-elevation populations would have only had to shift to nearby lower elevations via short distance dispersal to stay within their niche, contributing to greater stability and more tempered bottleneck signatures. By contrast, the northern range largely disappeared due to the presence of ice sheets (Shafer et al., 2010), alongside the most suitable connections to eastern B. vancouverensis populations (Lozier et al., 2013(Lozier et al., , 2016. Given that northern populations were likely much larger and well connected to other parts of the B. vancouverensis range prior to the LGM, as they are today, such habitat loss would have had disproportionately large effects on N e . Conversely, connections remade after the LGM could supplement N e via colonization from multiple founding populations (Ortego et al., 2015). Our results thus suggest that booms in heterozygosity from massive spatial and N e expansions in the north during favorable climates likely offset any busts from contractions in less favorable periods.
Focusing on bumble bees, a recent analysis of Bombus huntii in western North America  found nearly identical relationships for genetic diversity and Pleistocene niche stability, with northern populations exhibiting both relatively high contemporary heterozygosity and relatively low predicted LGM niche suitability compared with southern parts of the species range. Interglacial and postglacial increases in N e for species like B. vancouverensis and B. huntii contrasts with some expectations for cold-adapted bumble bees in other regions, which can expand their ranges during glacial periods and contract to interglacial refugia in warm periods (Hewitt, 2011;Martinet et al., 2018;Williams et al., 2018). However, models of reduced population size during glacial periods have also been supported for montane insects of western North America (Schoville & Roderick, 2009). Thus, diversity in bumble bees and other widespread montane species is likely determined by the degree of habitat stability alongside the arrangement of that habitat that can facilitate or constrain population expansion in particular regions in particular climates (Lee et al., 2019;Williams et al., 2018Williams et al., , 2022. For many cold-adapted North American Bombus species, southern parts of species ranges tend to be restricted to relatively high elevations. While such parts of the range may be relatively stable as populations can simply shift down and upslope to nearby areas with climate fluctuations (Lee et al., 2019), maximum habitable areas, and consequently Ne, may be reduced compared with high-latitude populations. Population dynamics should vary in parts of the globe with different orientations of major mountains, and it would be interesting to compare our results with other regions.
The potential connections of northern populations to other parts of the B. vancouverensis range discussed above raises an important caveat related to SMC models in general. Although employing multiple SMC methods can strengthen inferences, like any method results are sensitive to the assumed model (reviewed in Mather et al., 2020).
For example, periodic changes in migration and population structure like that discussed above can also influence the coalescent N e and the shape of SMC-inferred trajectories without major changes in true population size (Chikhi et al., 2018;Mazet et al., 2016). By comparing many individual genomes across a species range, consistent variation in the magnitude or timing of population size peaks and troughs still provide evidence for historical demographic changes influencing contemporary diversity. Because genomes from southern populations show persistently low N e compared with northern genomes, we infer that processes during the Quaternary established conditions that favor the maintenance of relatively low genetic diversity without needing to invoke strong effects of very recent (i.e., Holocene, Anthropocene) processes. The correlations with predicted habitable area from SDMs also provide support for the reliability of the SMC models; however, because we expect changes in population structure over time, some inferred differences could be due to migration rates rather than a simple reflection of the numbers of individuals alone. Another potential concern is that, while temporal SDMs have been commonly used for understanding patterns of genetic diversity, these models assume that a species niche is conserved through time, with no temporal adaptation to changing conditions that might allow persistence in regions with low predicted suitability based on contemporary occurrences. Local adaptation may be present across the B. vancouverensis range (Heraghty et al., 2022), and it is possible that both the SDM projection at the LGM is pessimistic and that N e trajectories could be altered by localized selection signatures in the genome, both of which could influence our conclusions. The temporal concordance of range size and N e dynamics provide some support for conclusions from each method, but the development of methods for coalescent inference and niche modeling that incorporate space, time, and natural selection would greatly assist our understanding of historical demography in B. vancouverensis and most other species (Chikhi et al., 2018).
We also cannot unambiguously state that recent demographic changes are not influencing B. vancouverensis alongside Quaternary fluctuations because SMC methods tend to perform relatively poorly at contemporary scales. However, the available San Juan Island genome  is a useful comparison to mainland samples by providing something of positive control for the speed with which SMC models can track diversity loss when isolation is more extreme. We expected that recent changes could be detectable in bumble bees, even for "long-term" SMC estimators, due to their annual generations. Indeed, unlike other samples, the San Juan Island genome indicated a striking contemporary N e decline. Although predicted suitability fluctuated over time and contemporary suitability is high, the island's isolation following the glacial retreat likely would have led to conditions favoring significant genetic isolation . Strong genetic signatures of isolation are common in bumble bees on coastal islands (Goulson et al., 2011;Jha, 2015;Lozier et al., 2011), providing a model of habitat fragmentation effects that may become common under climate change (Frankham, 1996). If lower elevational range limits erode under climate change (Kerr et al., 2015), populations in southern parts of the B. vancouverensis range will likely be pushed upslope and become increasingly isolated (Parmesan, 2006;Urban, 2018). This "sky island" isolation could begin to mimic the signatures observed on true islands (Epps et al., 2006;McDonald & Brown, 1992;Urban, 2018), even if long-term range dynamics drive most of the current variation in heterozygosity. Other works in bumble bees suggest Arctic populations may be at even greater risk than montane populations under climate change (Lee et al., 2019), so high-latitude B. vancouverensis populations are also potentially at risk despite their currently greater habitable area and genetic diversity. Additional work to understand demographic changes within the last century and monitor future genetic diversity will become important components of montane Bombus conservation.
In conclusion, demographic analyses suggest that populations of montane bumble bees may rapidly respond and maintain genetic diversity following climate-associated range perturbations and that gradients in diversity can be explained to a large extent by differences in these historical fluctuations. Some caution is thus likely warranted when attempting to link patterns of contemporary Bombus genetic diversity to human activities. To some extent, such patterns are encouraging for the long-term genetic health of bumble bee populations. However, there is reason to expect that future climate change will drive montane pollinator populations further upslope (Kerr et al., 2015), and given that long-term N e and heterozygosity appear constrained in high-elevation populations, warming temperatures may further erode habitable areas and population connectivity in this region. The distinct postglacial bottleneck inferred for the only truly insular sample suggests that population isolation is a concern for the erosion of genetic diversity when sufficiently strong and persistent. Thus, caution is also warranted for predicting future N e stability from relatively large N e ′s of the past. Additional sequencing of Bombus from other fragmented populations would provide insights into the temporal dynamics and consistency of reduced population sizes in insular habitats like islands, mountaintops, or disturbed habitat patches.

ACK N OWLED G M ENTS
We thank the National Science Foundation for support while conducting this research (URoL 1921585, DEB 1457645). We thank J Jackson, M Pimsler, and numerous undergraduate researchers for assistance with laboratory work. Permits used for the collection of

CO N FLI C T O F I NTE R E S T S TATE M E NT
Neither I nor my co-authors have a conflict of interest. No materials should require permission to reproduce.

DATA AVA I L A B I L I T Y S TAT E M E N T
All sequence data used in the study have been deposited at the National Center for Biotechnology Information Sequence Read Archive (see Appendix S1, Table S2 for accessions). Genetic data outputs, coordinate data from GBIF.org, filtered records and random background points with climate data, stability metric calculations, and example code are deposited on DRYAD (10.5061/dryad. dbrv15f4m). https://doi.org/10.5061/dryad.dbrv1 5f4m.